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ABSTRACT 

We present predictions drawn from cosniological hydrodynamic simulations for the 
physical, photometric and emission line properties of galaxies present during the latter 
stages of reionization from z = 9 6. We find significant numbers of galaxies that have 
stellar masses exceeding IO^Mq during this epoch, with metallicities exceeding one- 
thirtieth solar. Far from primeval "first-star" objects, these objects exhibit a significant 
Balmer break, are likely to have reionized their infall regions prior to z = 9, are 
dominated by atomic rather than molecular cooling, and are expected to be forming 
few if any metal-free stars. By z = 6, the space density of > 10^°Mq objects is 
roughly equivalent to that of luminous red galaxies today. Galaxies exhibit a slowly 
evolving comoving autocorrelation length from z = 9 — > 6, continuing a trend seen 
at lower redshifts in which the rapidly dropping bias counteracts rapidly increasing 
matter clustering. These sources can be marginally detected using current instruments, 
but modest increases in sensitivity or survey area would yield significantly increased 
samples. We compare to current observations of the z « 6 rest-UV and Lya line 
luminosity functions, and find good agreement. We also compare with the z ^ 7 
object observed by Egami et al., and find that such systems are ubiquitous in our 
simulations. The intrinsic Lya luminosity function evolves slowly from z = 9 — > 6, 
implying that it should also be possible to detect these objects with upcoming narrow 
band surveys such as DAzLE, if as we argue the detectability of Lya does not drop 
significantly to higher redshifts. We make predictions for near-IR surveys with JWST, 
and show that while a high density of sources will be found. Population III objects may 
remain elusive. We present and compare simulations with several recipes for superwind 
feedback, and show that while our broad conclusions are insensitive to this choice, a 
feedback model based on momentum-driven winds is mildly favored in comparisons 
with available data. 

Key words: galaxies: formation, galaxies: evolution, galaxies: high-redshift, cosmol- 
ogy: theory, methods: numerical 



1 INTRODUCTION 



The epoch of reionization is an important event in the his- 
tory of our universe, during which the light of the first stars 
and galaxies ate away at the neutral hydrogen fog and trans- 
formed the intergalactic medium (I GM) into the highly ion - 
ized cosmic web that we see today (iLoeb fc Barkana l200ll) . 
Detecting this epoch has long been a goal of observational 
cosmology. Recent observations of a Gunn-Peterson trough 
in z > 6 quasar spectra sug gest that the end of reioniza- 
tion epoch has now been seen iPan et al. ]^02). Conversely, 
third-year data from the Wilkin son Microwave Anisotropy 
Probe (WMAP;|Paie et al. suggests that the universe 

was predominantly ionized out to 2: ~ 10. When exactly 



reionization begins, how it proceeds, and what the sources 
of the ionizing photons are remain open questions. 

In light of this, a key goal of next generation ground 
and space-based telescopes is to detect the objects re- 
sponsible for reionizing the universe at z > 6. At 
present these are believed to be star-forming galaxies since 
there are insuffi cient quasars at z > 6 to do the job 
jPan et al. 11200 ih . although the latest SDSS samples show 
a faint-end slope that is tan taUzingly close to the required 
steepness JPan et al. I [200^. Pew reionization-epo ch ob- 
jects ha ve been observed up till now Ij lhoads & Malhotr^ 
1*2004; Ega mi et al"1l2005l: IStern at al. 120 05: Hu et al. 200^; 
ICharv. Stern fc Eisenhardt Il2005^ ■ and hence their detailed 
physical properties remain poorly constrained. In order to 
optimize detection strategies with current and future facili- 
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ties, it is important to understand the nature of these early 
galaxies. 

Simulations of galaxy formation can provide insights 
into the expected properties of these systems within the 
context of a now well-established framework for hierarchi- 
cal structure formation. While the properties of dark matter 
halo assembly are relatively well understood owing to pre- 
cisely dete rmined cosmological parameters (ISpergel et al. I 
l2003l . libofit) . the properties of the gas clouds that condense 
to form observable stars remain subject to a large range of 
theoretical uncertainties. Previous studies have utilized an- 
alytic or semi-analytic approaches to predict the properties 
of early galaxies, based o n simplified models for the rele- 
vant baryonic physics (e.g. iBarkana fc Loeb1l2nnil and ref- 
erences therein). An alternative approach is to employ full 
cosmological hydrodynamic simulations that account for all 
the relevant cooling and star formation processes associated 
with galaxy formation, as presently understood, constrained 
by rapidly-improving observations of high-redshift systems. 
Such simulations carry information about the full three- 
dimensional structure of accretion and feedback processes 
that are vital to regulating galaxy formation, and track the 
growth of baryonic galaxies within a full hierarchical struc- 
ture formation scenario. These simulations have been fairly 
successf ul at modeling post-reionization high-redsh ift galax- 
ies (e.g. iNight et al.T 120051: iFinlator et aL]l2005^ ■ so they 
provide a natural framework to extend studies back into 
the reionization epoch. Such an a pproach is complementary 
to mod els of first s tar formation llAbel. Brvan. fc NormarTI 
[2000 ; B romm fc La rson 20o3) at 2 > 15, which aim to pre- 
dict the properties of primeval galaxies by understanding 
early stellar assembly from primordial gas. 

A significant difficulty with using models or simula- 
tions tuned to match lower redshifts is that the physics of 
galaxy formation may be quite different during the epoch 
of reionization. Compared to today, the universe during 
the reionization epoch had higher densities, shorter cool- 
ing times, lower metallicities, and a fiatter power spec- 
trum on mildly nonlinear scales, all of which could in- 
fluence the formation and early evolution of galaxies. 
Furthermore, the process of reionizatio n itself may have 
large effects on the galaxy population iBarkana fc Loebl 
120011) . perhaps supp ressing star formation through pho- 
toionization heating jB arkana fc Loeb ' 'l999*) or stripping 
jSigward. Ferrara. fc Scannapieco .2005.) , or perhaps en- 
hancing it through shock-induced star formation or through 
early black holes emit ting X-ra ys that stimulate molecular 
hydrogen production jHaimanfe oSl. It is therefore possi- 
ble that our recipes and models for galaxy formation tuned 
to observations at the present epoch may have little rele- 
vance during the reionization epoch. Indeed, this is almost 
certainly the case for primeval galaxies, i.e. the very first 
galaxies forming nearly metal-free stars. 

On the other hand, if we want to understand the prop- 
erties of galaxies that will be observable in the near future 
with next-generation instruments such as JWST and 20- 
30m class telescopes, it is questionable whether these first 
galaxies forming their very first stars are the relevant ob- 
jects to investigate. In this paper we argue that galaxies 
that will be detected in the foreseeable future have already 
formed many generations of stars, as they are among the 
largest and/or most vigorously star forming galaxies at these 



epochs. They have circular velocities well above that ex- 
pected to be strongly affected by photoionization, metallic- 
ities well above the threshold for a "normal" stellar initial 
mass function, and halo virial temperatures well above the 
threshold where atomic lines become the dominant coolant. 
Furthermore, these systems are highly biased, which means 
that they are expected to be relatively unaffected by the 
global reionization process occuring in the universe at large, 
because they likely reionized their local volume at an earlier 
epoch. Hence we argue that observable galaxies at 2 < 9 are 
probably more similar to lower- redshift systems than they 
are to early primeval systems, and that simulation methods 
developed to model post-reionization galaxy formation may 
provide a reasonably realistic description for these systems 
that are responsible for completing the process of reioniza- 
tion. 

In this paper we study the detailed properties of the 
galaxy population at 2: ~ 6 — 9 using cosmological hydro- 
dynamic simulations. We include physical processes such 
as kinetic feedback with accompanying metal injection and 
metal-line cooling, in addition to the usual radiative and star 
formation processes required for galaxy formation, while ig- 
noring the highly complicating effects of radiative transfer. 
Our primary goal here is to use such models to make pre- 
dictions for reionization-epoch galaxies that can be tested 
against present and upcoming observations. We find that 
many galaxies with stellar masses exceeding 1O*M0 are al- 
ready in place at 2 = 9, and by 2 = 7 there are plenty with 
M* > 10^ Mq. These galaxies have all been enriched at sig- 
nificant levels, to at least 3% solar, and more typically 10% 
solar. They are highly biased and highly clustered, and sim- 
ple estimates show that they will reionize their infall regions 
by 2 ~ 9 even if the rest of the universe remains mostly 
neutral. They would be detectable in large numbers with 
near-infrared surveys only modestly deeper than accessible 
today, and depending on line transfer effects could also be 
detectable in Lya emission, perhaps even more abundantly 
than with broad-band surveys. We compare with observed 
2 ~ 6 — 7 systems, and demonstrate that our models nicely 
reproduce available observations. In short, the aim of this 
paper is to set the stage for detailed comparisons of 2 > 6 
objects with predictions from hierarchical models of struc- 
ture formation including gas physics, and show that initial 
comparisons can already provide valuable insights into the 
nature of galaxies in the reionization epoch. 

Our paper is organized as follows: In ij^jwe present our 
simulation methodology, details of our various runs, and al- 
gorithms for determining galaxy properties. In SjSlwe present 
physical properties of 2 = 6 — 9 systems, including their 
masses, star formation rates, metallicities, and halo prop- 
erties. In S0] we discuss galaxy clustering and its evolu- 
tion, along with implications for local reionization. In JS] 
we present predictions for broad-band and Lya emission 
properties of high-redshift galaxies, including comparisons 
to available observations and prospects for detection with 
future instruments. We summarize our results in ^ 



2 SIMULATIONS 

We employ the parallel cos mological galaxy formation code 
Gadget-2 iSpringel Il2005^ in this study. This code uses an 
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''Additionally, a suffix "nw" , "cw" , "mzw" is added to 
denote no, constant, or momentum winds, respectively. 

''Box length of cubic volume, in comoving /i~^Mpc. 

'^Equivalent Plummer gravitational softening length, 
in comoving h~^kpc. 

"^AU masses quoted in units of 10^ Mq. 

''Minimum resolved galaxy stellar mass. 

entropy-conservative formulation of smoothed particle hy- 
drodynamics (SPH) along with a tree-particle-mesh algo- 
rithm for handling gravity. It includes the effects of radia- 
tive cooling assuming ionization equilibrium with primor- 
dial chemistry, which we have extended to include metal-line 
c ooling based on the collisional ionization equilibrium tables 
of ISutherland fc Dopita I il993l) . The metal cooling function 
is interpolated to each gas particle's metallicity as tracked 
self-consistently by Gadget-2, and added to the primordial 
cooling rate. Gas particles that are eligible for star forma- 
tion are continually enriched based on an instantaneous re- 
cycling approximation with a yield parameter of 0.02. The 
thermal energy of each gas particle is evolved on its cool- 
ing timescale (if it is shorter than the Courant or dynami- 
cal timcscalo) assuming isochoric conditions. For full details 
sec Qppcnheim cr & Dave (200(i). Star s are f ormed using a 
recipe that reproduces the iKennicutt I ^9^ relation, em- 
ploying a subresolution multi-phase model that tracks con- 
den sation and evaporati o n in th e inte rstellar medium follow- 
ing ^McKcc fc Ostriker I (^T^); see 'Spri ngel fc Hernguist I 
i2003 a) for details. Stars inherit the metallicity of the par- 
ent gas particle when spawned, and from then on cannot be 
further enriched. 

Gadget-2 has an observationally-motivated prescrip- 
tion for driving superwinds out of star forming galaxies, 
which is tuned to repro duce the global stellar mass d en- 
sity at the present epoch iSpringel fc Hernguist ll2003tJl . In 
lOppenheimer fc Dave I tOOdi) we describe modifications to 
this scheme designed to better match observations of IGM 
enrichment. The free parameters are the wind velocity of 
material ejected and the mass loading factor, which is the 
ratio of the rate of matter expelled to the star formation 
rate. Here, we will concentrate on three superwind schemes. 
The first is one where superwinds are turned off (though 
thermal feedback from supernova is included). This "no 
winds" (nw) model is known to ov erproduce the amount 
of star s significantly at low redshifts iSpringel fc Hernguist I 
lioolB), so its purpose here is to gauge the effects of the 
other (more observationally consistent) superwind models. 
The second is a "constant wind" (cw) model where all the 
particles entering into winds are expelled at 484 kms~^ 
out of star forming regions and a constant mass loading 
factor of 2 is assumed, which is the scheme used in the 



runs of ISprineel fc Hernguist 1 ll2003d) . Finally, we intro- 
duce a new "momentum wind" (mzw) model where the 
wind velocity is proportional to the local circular veloc- 
ity («m = 3 * Vcivc, where Ucirc is computed from the lo- 
cal gravitational potential), and the mass loading factor 
is inversely proportional to the velocity dispersion (77 = 
300 km s~^/g), as expected in momentum-dri ven wind sce- 
narios iMurrav. Quatert. fc T hompson 2005) and consis- 
tent with observations bv iMartin (_20P5J. The mometum 
wind model also gives a velocity boost in low-metallicity 
systems, based on the arguments that more UV photons 
are produced per unit stellar mass at lower metal licities 
(specifically, we employ eqn. 1 of ISchaerer I l2003h . and 
that it is these UV photons that ar e driving the wind 
jMurrav. Quatert. fc Thompsoni |200^. Note that these 
supernova-driven winds must be included "by hand" mainly 
because we lack the numerical resolution to self-consistently 
follow the physical processes required to drive bulk flows out 
of star forming regions. 

These superwind models are primarily intended to il- 
lustrate the sensitivity of our results to our somewhat ad 
hoc feedback prescriptions, though admittedly one could en- 
vision a wider range o f feedb ack models. We note that in 
lOppenheimer fc Dave I ll200d) we find the momentum wind 
model to be in better agreement with a range of observations 
including intergalactic C iv absorption from 2 ~ 5 — > 1.5, 
and our comparisons in this paper also mildly favor the mo- 
mentum wind model. We further note that observed typical 
outfiow velocities for z ~ 2 — 3 Lyman break galaxies of 
hundreds of km/s JPettini et al. 12 001: Sliaplcv ot al. 200j) 
and t heir inferred mass loading factors of r; 4 IIErb et al. I 
l2006f) are both consistent with our momentum wind model 
parameterization. For these reasons we generally focus on 
the momentum wind model when the conclusions from the 
different wind models are broadly similar. 

Gadget-2 implements heating due to a photoionizing 
background, which we t ake from the latest CUBA model 
jHaardt fc Madau llioOlh that assumes a 10% escape frac- 
tion of ionizing radiation from galaxies in addition to the 
contribution from quasars. For this model, reionization oc- 
curs at 2 « 9. Spatial uniformity is assumed, which is obvi- 
ously incorrect in detail during reionization, but may nev- 
ertheless be a passable local approximation if the galaxies 
of interest reionize their surroundings earlier and if their 
masses are significantly above that which would be sup- 
pressed by photoionization. For galaxies resolved by our 
simulations, we will argue that both those statements are 
typically true. 

We run nine simulations: Three different volumes as 
listed in Table each with our three different superwind 
schemes. These three volumes can be used to judge numer- 
ical resolution effects at overlapping mass scales. We note 
that Gadget-2 has been extensively tested for resolution 
convergenc e with regards to its star formation and feedback 
algorithm iSpringel fc Hernguist Il2003at iNight et aL I I2OO5I: 
iFinlator et al. Il2005^ ■ so we do not conduct a wide suite of 
resolution tests. Nevertheless we are reassured by the good 
numerical convergence seen among our runs here, as we will 
point out in upcoming sections. The choice of box size has 
also been carefully considered, because overly small volumes 
can lead to galaxy format ion occuring too late owi ng to miss- 
ing power on large scales (iBarkana fc Loeb l2004l) . while too 
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large volumes would be unable to resolve the early objects 
of interest. Our box sizes are therefore chosen in a range 
spanning between these pitfalls. 

All runs use cosmological parameters consistent with 
recent joint analyses of first-year WMAP, S PSS, Lyman al- 
pha f orest, and Type la supernovae data ijTeemark et al. | 
I2OO4I: ISeliak et al.ll2n04h : We assume = 0.3, A = 0.7, 
Ho = 70 kms"^ Mpc"\ erg = 0.9, and = 0.04. More 
recent third-year WMAP data indicates a lo wer amplitude 
of th e matter power spectrum {as ~ 0.75; ISoergel et al. I 
|200(J) . which would result in fewer galaxies, but also a lower 
Q and higher ^6 which would result in more galaxies at 
these epochs . Howe ver, there is still some debate on as, as 
ISeliak et al. I (l200dl determines a significantly higher value 
(o"8 = 0.85) by incorporating Ly-a forest constraints. We 
expect that our results would therefore be qualitatively sim- 
ilar with the correct cosmology, and indeed there are signifi- 
cantly larger uncertainties associated with modeling details 
such as our superwind scheme. 

Each of our runs heis 256^ dark matter and 256"^ gas 
particles, evolved from the linear regime down to z — 6. 
The smallest volume is a cube of 8/i~^Mpc on a side, with 
an equivalent Plummer gravitational softening length of 
0.625/i~^kpc (comoving; i.e. better than 100 pc physical at 
2: > 6), and our largest volume has a 32/i~^Mpc box length. 
The initial conditions are chosen t o be random subvolumes 
of the universe, generated using an iEisenstein &: Hu I lll999ll 
power spectrum when the universe is still well within the 
linear regime (typically 2 > 150). Since the nonlinear scale 
at 2 = 6 is less than our smallest box length, each volume 
should contain a representative galaxy population, though of 
course the rarest objects will be absent. Given the canonical 
value of 60 particles for robu st identification of dark matter 
halos jWeinberg et al.^ll999^ ■ our smallest volume is resolv- 
ing halos with mass > 2 x 10* M0, with an associated virial 
temperature > 10* K. As we will show in Hlj.5l halos con- 
taining resolved galaxies (defined below) are substantially 
larger, typically above 3 x 10^ M©. This means that primor- 
dial atomic gas should be the dominant coolant in all our 
resolved systems. Note that our simulations do not include 
molecular hydrogen cooling. 

We identify galaxies using Spline Kernel Interpolative 
DENMAX, and dark matter halos using t he spherical over- 
density algorithm (see iKeres et al. I |2p05!, for full descrip- 
tions). We only consider galaxie s with stellar masses ex- 
ceeding 64 star particles, which in lFinlator et al.H ((200^ we 
determined represents a converged sample both in terms of 
stellar mass and star formation history. The minimum re- 
solved stellar mass of galaxies in each volume is listed in 
Table In order to fully utilize the dynamic range in our 
three simulation volumes, we present stellar mass functions 
derived by splicing together the results in our various vol- 
umes, in the following manner: For masses where multiple 
volumes satisfy the 64 star particle criterion, we take the 
maximum value among our simulated mass functions. In 
practice, at any mass this is usually the mass function of 
the largest volume (or close to it), and hence we effectively 
present mass functions from the largest volume that resolves 
a given galaxy mass. 

The resulting 2 = 6 stellar mass functions are shown 
in Figure for the momentum wind (top panel), constant 
wind (middle) and no wind (bottom) models. The dotted. 





1 1 1 1 1 1 1 1 1 1 


mzw model _ 
2.2,0.2% : 

_ 


-4 


1 1 1 1 1 1 1 1 1 








1 1 1 1 1 1 1 
cw model ! 
3.0,0.6% _ 




1 '^r 


-4 


1 1 [ 1 1 1 1 1 1 








1 1 1 1 1 1 1 
nw model \ 

4.8,1.3% : 


-4 


^ 1 1 - 1 1 1 


1 1 1 1 1 1 1 1 1 



7 8 9 10 11 

log M. (MJ 

Figure 1. Stellar mass functions of galaxies at 2 = 6 in the 
momentum wind (mwz) model (top panel), constant wind (cw) 
model (middle panel), and no wind (nw) model (bottom panel). 
Each panel shows the aggregate mass function (solid line) as as- 
sembled from the individual mass functions in the 8/i~^Mpc box 
(dot-dashed line), 16/i~^Mpc box (dashed line) and 32h~^Mpc 
box (dotted line). Vertical lines of the same line type indicate the 
64 star particle mass threshold for each volume. The numbers in 
the upper right represent the percentage of baryons in galaxies 
and in stars, respectively. The mass functions join smoothly onto 
each other, showing good resolution convergence in the overlap 
region, and lending confidence to the procedure of stitching to- 
gether the mass functions. 



dashed, and dot-dashed curves represent the mass functions 
from the 32, 16, and 8/i~^Mpc volumes. Overlaid is the ag- 
gregate mass function, shown as the solid line, which as ex- 
pected tracks the 32/i~'^Mpc mass function at high masses 
and the 8/i~^Mpc mass function at small masses. The ver- 
tical lines represent the 64 star particle mass limit for the 
three volumes, which is where the turndown begins, at least 
for the feedback models; for the no wind model the turndown 
begins at slightly larger masses. The caption in the upper 
right shows the mass fraction of baryons in cold gas and 
stars at 2 = 6 in these three models, showing that the mo- 
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mentum wind model produces significantly less cold gas and 
stars than the constant wind model by this epoch, which in 
turn produces significantly less than the no-wind case. The 
key point from Figure is that mass functions of the differ- 
ent resolution simulations agree quite well in the overlapping 
mass range, resulting in a nearly seamless mass function ex- 
tending from lO"^-^ Mq to beyond 1O^°M0. Hence we do not 
expect that numerical resolution limitations will affect our 
results, so long as we concentrate on galaxies exceeding our 
64 star particle mass criterion. 

To obtain photometric properties of our simulated 
galaxies, we use [Bruzual & Chariot (2003) population syn- 
thesis models to convert the simulated star formation history 
of each galaxy into broad band colours in an assortment of 
filters, including J, Ks, and the Spitzer IKAC [3.6/i], [4.5/i], 
[8fi] and MIPS [24/i] bands (we use AB magnitudes through- 
out). We describe this further in H5.'2\ We account for dust 
reddening using a prescription based on the galaxy's metal- 
licity and a locally-calibrated metallicity-extinction relation. 
We accou nt for IGM abs orption bluewards of rest-frame Lyo 
using the lMadau I lll995t) prescription; at these redshifts this 
results in virtually complete absorption, so we do not con- 
sider results fro m bands at wavelengt hs shorter than red- 
shifted Lya. See lFinlator et cil~l (|200a) for a more complete 
description of these procedures. We also predict Lya emis- 
sion line properties based on the instantaneous star forma- 
tion rates in each galaxy as output by Gadget-2; we de- 
scribe this in more detail in il5.8l 

We concentrate on studying galaxies between z = 6 ~9 
for several reasons. First, our simulations have limited dy- 
namic range, which means that pushing to earlier epochs re- 
sults in a statistically insignificant sample of resolved galax- 
ies that have formed by that epoch (particularly in our larger 
volumes). Second, our assumed ionizing background turns 
on at 2; ~ 9, which means that prior to that time the sup- 
pression of gas infall due to local photoionization, which may 
be relevant, is not taken into account. Finally, there is the 
practical consideration that surveys of reionization galaxies 
are likely to concentrate on galaxies in this epoch because 
they are both nearer and intrinsically brighter than those at 
higher redshifts. 



3 PHYSICAL PROPERTIES 

3.1 Global Star Formation Rate Density 

We begin by showing some basic physical properties of 
reionization-epoch galaxies in our simulations. The global 
star formatio n rate density evolu tion, sometimes called the 
Madau plot iMadau et al. llT996h . has become a bellwether 
test for models of cosmological galaxy evolution. Figure |5| 
shows the global star formation rate density for the three 
wind models. These values were computed from the resolved 
galaxy population in the three vol umes (8,16,32/t~^Mpc) 
for eac h wind model. As found by ISpringel fc Hernauist I 
(lioosB), superwind feedback has a large effect on sup- 
pressing star formation in the universe, even at these 
early epochs. The constant wind case shows increasing 
suppression with time relative to the no wind scenario, 
while the momentum wind model provides a relatively 
redshift-independent global suppression of ~ x 10. Note that 








O 



Redshift 

Figure 2. Global star formation rate density (Madau plot) for 
the momentum wind (mzw; solid line), constant wind (cw; dotted 
line), and no wind (nw; dashed line) models. Open data point s 
with error bars are taken from a compilation bv lHopkins I J200^, 
which includes a star formation rate-dependent extinction cor- 
rection. Open squares are other data point s at z: ^ 5 .5 from the 
literature: _Fontana ct al. 1 20 oj, z = 5.5), iBunker e t al. 1 20oJ, 
low point at 2 6), |Bouwe ns et al. 1 120061 . high point at 2 f» 6), 
and iBouwens et al. rl20o3 " 2 ^ 7.S) . The constant wind model 
causes a suppression of fsi x3 versus no winds by z = 6. The mo- 
mentum wind model provides further suppression by ~ x3 — 4 
(depending on redshift) versus the constant wind case. Extrapola- 
tions of both constant and momentum wind models appear to be 
in agreement with observations, while if no superwind feedback 
is implemented there is significant overcooling. 



our simulations includ e metal -line cooling, unlike those of 
ISpringel fc Hernaxiist I (|2003lJ) . but at these early epochs 



primordial cooling is so efficient that the additional cool- 
ing provided by metal lines does not significantly impact 
the star formation efficiency jHernauist fc Springel I l2003l : 
iKeres et al. "20051). 

Fi gure 121 also shows data points from a compilation 
bv llfopkinsi y004) (only the z > 3 data are shown). We 
have supplemented this with more recent measurements at 
2 > 5.5 as noted in the caption. Comparing the simula- 
tions to the data, it is evident that superwind feedback is 
required in order to bring the simulations into broad agree- 
ment with observations and s olve the overcooling problem 
JSprineel fc Hernauist 1l2003bl) . Based solely on these data 
it is not clear which of our superwind models produces best 
agreement; indeed a simple extrapolation of our two models 
seems to roughly bracket the observed range at z ~ 3 — 6. 
The exceptions are the z ~ 6 p oint from [ Bun ker et alH 
(120041!) ■ and the z « 8 point from iBouwens et al. I J2004l) . 
which both lie significantly below the predictions. However, 
given that these points are based only on a handful of bright 
objects detected at those redshifts, it is unclear what to 
make of this discrepancy. In general, observational selection 
effects could significantly impact this comparison, as the 
galaxies identified in our simulations may not be equivalent 
to those used to determine the data points. We leave more 



6 Dave, Finlator, Oppenheimer 







1 1 


1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 

mzw model _ 










z = 6 


- 




n 






z = 7 










z = 8 


- 










z = 9 


- 


<J — 


2 












A 








- 

N — 






-1— |— 


1 1 1 1 1 1 1 


— 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 


^ — 1 — 1 — 1 — 

^ ~ 
Z = D _ 










mzw 


- 









~ - 


cw 










~~r--~^^ nw 


- 














<j — 


























4 








\ - 






-1— |— 


1 1 1 1 1 1 


1 1 1 1 1 1 


1 1 1 1 
z = 6 ! 










mzw 













cw 










^"^^^ ^ nw 
















O- 


2 




















X 

\ 

\ 












■ . \ 




4 


~ 1 


barypmc 


1 


" \ - 



J I I I LJJ I I I I I I I I I I I 1 I I I 

7 8 9 10 11 

log M (MJ 



Figure 3. Top; Stellar mass functions of galaxies for ^ = 9, 8, 7, 6 
from the momentum wind model. Units are number per cubic co- 
moving ?t~^Mpc per log stellar mass. Vertical line is our galaxy 
stellar mass resolution limit of 64 star particles. As expected the 
stellar mass function increases with time, and it additionally ob- 
tains a bend such that the faint-end slope is somewhat shallower 
with time. Middle: Stellar mass functions for the mzw (solid line), 
cw (dotted line), and nw (dashed line) models at 2 = 6. This 
shows that the momentum wind model is most effective at sup- 
pressing star formation across all galaxy masses. Bottom: Same 
as middle panel, for total cold baryonic mass function (i.e. star 
forming gas -|- stars). The three wind models have much more 
similar baryonic mass functions as compared to their stellar mass 
functions. 

detailed comparisons for future work. The main point here 
is that superwind feedback is a necessary ingredient for sup- 
pressing early galaxy formation to obtain broad agreement 
with observed cosmic star formation rate density evolution. 

3.2 Mass functions 

Hierarchical structure formation models generally predict a 
rapid early buildup of stellar mass in galaxies. Stellar mass 
functions at 2: = 9 — » 6 are shown in the top panel of Fig- 



ure 13 for the momentum wind model. As expected, this 
function grows with time, such that there are about twenty 
times more galaxies with M ~ 10^ at z = 6 as compared 
to 2 = 9. The shape is initially close to a power law, and ob- 
tains a more pronounced bend with time; the faint-end slope 
of a Schecter function fit goes from —2.17 at « = 9 to —1.85 
at z = 6. A significant number of galaxies with stellar masses 
exceeding 10* Mq, and a few exceeding 10® Mq, are already 
in place by z = 9. By z = 6, galaxies exceeding 10^° have 
a space density of about 0.0002 per cubic /i~^Mpc. This is 
comparable to t he space density of t he SDSS Luminous Red 
Galaxy sample iZehavi et al. ll2'005^ . suggesting that these 
objects at 2: = 6 are the progenitors of large ellipticals at 
low redshift. 

The middle and bottom panels show a comparison of 
stellar mass functions and total galaxy baryonic mass func- 
tions for the various wind models. The superwind feedback 
prescription has a significant effect on the stellar mass func- 
tion, producing nearly a factor of three reduction in the mass 
function for the momentum wind model as compared with 
the constant wind model. The suppression of star formation 
in feedback models relative to the no wind case is actually 
larger at larger masses, counterintuitive to t he notion that 
feedb ack preferentially affects small galaxies (|Dekel fc Silk I 
Il986h : this reflects the hierarchichal assembly of large galax- 
ies from smaller ones. There is much less of a difference 
among models for the total baryonic mass function, because 
the constant wind model has a relatively low mass loading 
factor that does not expel much of the gas, while the mo- 
mentum wind model has lower wind velocities that allow 
fairly quick re-accretion. The differences seen here can be 
regarded as an estimate on the modeling uncertainty from 
our feedback algorithms. 

3.3 Stellar Birth Rates 

Figure 0] shows the star formation rates as a function of 
stellar mass at 2: = 6 and 2 = 8 in the momentum wind 
model. The three clumps of points correspond to galaxies 
from our three volumes, which smoothly connect with each 
other. Star formation rates are essentially proportional to 
the stellar mass at both redshifts, with the correlation be- 
ing marginally tighter at 2 = 8, implying a fairly constant 
stellar birth rate across all masses. This tight relationship 
between stellar m ass and star formati on rate is also seen 
at lower redshifts iFinlator et al. ll2005^ ■ and arises because 
star formation in high-redshift ga laxies is dominated by gas 
infall on dynamical timescales iKeres et al. I l20o3) rather 
than more stochastic merger-driven burs ts as in some semi- 
analytic models (e.g. iKolatt et ai~lll999l) . Given that these 
simulations cannot resolve the detailed internal dynamics 
of mergers that lead to dramatic short-lived starbursts, it 
is possible that the scatter in this relationship is underes- 
timated. Nevertheless the basic trend that the most rapid 
star formers are generally the largest galaxies remains a fun- 
damental prediction of hydrodynamic simulations of galaxy 
formation. 

Galaxies at a given stellar mass are actually forming 
stars approximately twice as fast at 2 = 8 than at 2 = 6. 
Note that the total star formation rate in resolved galaxies 
increases by a factor of three from 2 = 8 to 2 = 6 since there 
are many more resolved galaxies at 2 = 6. These statements 
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Figure 4. Star formation rate as a function of stellar mass in 
tile momentum wind model at z = 8 (open circles) and z = 6 
(crosses). The star formation rate and stellar mass are generally 
proportional; the dotted line shows a linear relationship for com- 
parison. 
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are also generally true for the constant and no wind models 
(not shown). The largest galaxies at both redshifts are form- 
ing stars in excess of 10 Mq/jx, which is vigorous even by 
today's standards. This indicates that the epoch of galaxy 
growth is already well underway by z ~ 8. 

The tight positive correlation between stellar mass and 
star formation rate is a basic prediction of hydrodynamic 
simulations at high redshifts that has yet to be verified. Ob- 
servations by iPapovich et al. I |20oF) suggest that at red- 
shift z ^ 2, such a relation is broadly in place, although the 
scatter is much larger than observed in simulations. Such 
a trend is certainly not observed today, as massive galaxies 
presently are forming stars very slowly if at all. It remains 
to be seen if simulations can quantitatively reproduce this 
so-called galaxy downsizing, though the qualitative trend is 
a n atural outcome in hierarchical galaxy formation models 
(see iDave et al. Il200l for more discussion). 

The importance of mergers in driving galaxy prop- 
erties at high redshifts remains hotly debated. In addi- 
tion to the issue of stochasticity in star formation due 
to mergers as discussed above, there is also the issue of 
feedback due to merger-related processes, including the 
impac t of central black hole growth and asso ciated feed- 
back llDi Matteo. Soringel. fc Hernauisnl2005^ . AGN feed- 
back has the potential t o drive out large quantities of gas 
from g alactic potentials dSpringel. Di Matteo. fc Hernauisti 
l2005lj) . an effect we are are not modeling here. It is plau- 
sible that at these high redsh ifts, where the quasar popula- 
tion is small (iFan et al. the impact of such processes 
will be minimal. However, this bears future testing, and is 
noteworthy as another caveat in modeling the evolution of 
reionization-epoch galaxies. 
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Figure 5. Top: Gas-phase metallicity versus stellar mass at z = 8 
(circles) and 2 = 6 (crosses) in the momentum wind model. The 
three "clumps" of points correspond to the three different sim- 
ulation volumes. There is little evolution in the mass- metallicity 
relation at these early redshifts. Line s hows the median mass - 
metallicity relation at z ~ 0.1 from SDSS iTremonti et al. 
Bottom: Gas-phase metallicity of galaxies at z = 8 versus stel- 
lar mass, for the momentum wind (circles) and constant wind 
(crosses) models. The factor of three difference is due to a similar 
difference in the efficiency of converting gas into stars, as can be 
seen in the bottom panel of Figure l3l 



3.4 Metallicities 

Metal pollution has a large effect on the reionization epoch, 
both in terms of the cooling curve th at strongly affects 
the Jeans mass in star forming regions l|Bromm fc Larson I 
^004), as well as enabling the formation of dust that can 
affect the detectability of systems, particularly in the Lya 
emission line. In the top panel of Figurel^we show the mass- 
averaged gas-phase metallicity of galaxies in the momen- 
tum wind model at 2; = 8 and 2 = 6. The mass-metallicity 
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Figure 6. Percentage of stellar mass formed in our 
simulations with metallicity below 10~* solar, as a function of 
time. The fraction is always quite low, well below 1% at all epochs 
in the wind simulations, and the fraction generally falls with time. 
This suggests that Population III star formation will be rare at 
2 ~ 6 — 9, and even rarer at later epochs. 



relation shows little redshift evolution during this epoch; 
as galaxies form more stars, they commensurately increase 
their metallicity so as to preserve a linear relation. For com- 
parison, the SPSS median mass-metallici ty relation is shown 
as the solid line jTremonti et al. Il2004() : interestingly, the 
slope is similar, but the amplitude is almost 1 dex lower 
than observed at z ~ 0. This means that although little 
evolution is detected from z — 8 6, the metallicity at a 
given stellar mass must increase by ~ x5 — 10 by redshift 
zero. Given the protracted period of time from z = 6 — > 0, 
this does not seem unre asonable, but still bears future exam- 
ination. More recently, lErb" et al. I i2006() has measured the 
meiss-metallicity relationship at z « 2, and found it to lie 
only 0.3 dex below the local relation; this is consistent with 
the relatively slow evolution predicted in these simulations. 

The bottom panel of Figure|S]shows the gas metallicities 
in the two wind models at 2 = 8. The feedback mechanism 
has a significant effect on the metallicity, with the mome- 
tum wind model having typically half a dex lower metallicity 
than the constant wind case. This is directly a reflection of 
the ~ X 3 lower amount of stellar mass in galaxies of a similar 
total baryonic mass, as shown in Figure^] While we broadly 
favor the momentum wind model, we note that the factor 
of 3 variation in the mass-metallicity relation amplitude is 
probably indicative of current feedback modeling uncertain- 
ties. We have not shown the stellar metallicities, but they 
are similar to the gas metallicities, being lower by at most 
a factor of two. 

The main point of Figure |5] is that in all cases, the 
metallicities of resolved galaxies exceed about one-thirtieth 
solar. While this would be considered metal-poor by to- 
day's standards, it still far exceeds the metallicity thresh- 
old of ~ 10"* that is purported to enact the transition to 
a normal Population II stellar population from a "metal- 
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Figure 7. Fraction of halo mass in total baryons (upper band) 
and stars (lower band) as a function of halo mass at z = 6 in the 
momentum wind (top panel), constant wind (middle panel) and 
no wind (bottom panel) models. With no outflows, halos contain 
more than the global fraction of baryons, indicated by the hori- 
zontal dashed line. Outflows reduce the baryon fraction relative 
to the no wind case by around 30% and 20% in the momentum 
and constant wind models, respectively. The stellar fractions a 
stronger trend with halo mass, reflecting primarily the later col- 
lapse time of smaller halos. The momentum wind case shows a 
signiflcant suppression of star formation relative to the other two 
models, owing to the high mass loading factors in these small 
early systems. 



free" Population III one (see iBromm fc Larsorni2004j . and 

references therein) . Hence it is expected that these high red- 
shift galaxies should already have little if any contribution 
from Population III stars, assuming that the metal mixing 
timescale in the ISM is shorter than a few dynamical times. 
Furthermore, these galaxies may also contain a significant 
amount of dust if the timescale for dust formation is as short, 
as has been inferred f rom high-redshift quasar observations 
jMaolino et al. Il2004) . 

While the mean metallicity in galaxies is large compared 
to the canonical value required to suppress Population III 
star formation, it is still possible that some stars continue to 
form with very low metalli cities at late times. For instance, 
Ijimenez fc Haiman I (l2006l) interpreted recent data from z ^ 
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3 — 4 galaxies as implying significant (10 — 30%) ongoing 
Population III star formation. If verified, this would suggest 
that z > 6 galaxies should be forming even greater fractions 
of low-metallicity stars. 

In Figure |S| we calculate the percentage of stellar mass 
formed in our 8/i~^Mpc simulations having Z < W~^Zq as 
a function of time. We choose our smallest volume to maxi- 
mize the likelihood of early, metal-free stars; the values from 
larger volumes were found to be lower. The percentages are 
nonzero but very low, typically less than 0.5% at all epochs, 
with the only exception being the no wind case where un- 
fettered infall leads to an initial burst of near-metal free 
star formation (still only ~ 1.5%). The stochastic nature of 
our star formation algorithm, and the fact that we are not 
forming individual stars but rather star "clumps" that are 
means that the values shown should be re- 
garded as crude estimates. However, it is still clear that our 
simulations produce very little Population III star forma- 
tion at 2 = 6 — 9, likely dimishing furth er to lower redshifts, 
and ar e therefore inconsistent with the I Jimenez fc Haiman I 
i2006f ) hypothesis. Future observations will be required to 
determine if this represents a failing of our models. 

3.5 Halo Baryon Fractions 

During the reionization epoch, galaxies are likely embed- 
ded in a inhomogenously percolating radiation field. Due to 
the computational complexity of radiative transfer, our sim- 
ulations do not track this, and instead assume a uniform 
photoionizing background turned on at 2 ~ 9. Clearly this 
approximation is incorrect in detail; however, we can as- 
sess whether the resulting galaxy population is likely to be 
significantly affected by the external ionizing radiation by 
examining some basic properties of baryons in halos. 

Figure |7| (top panel) shows the mass fraction of total 
baryons (upper band of points) and stars (lower band) in ha- 
los that host galaxies above our 64-particle stellar mass res- 
olution limit, for our momentum wind runs. There are three 
clumps of points corresponding to halos in our three vol- 
umes; the agreement between the runs at overlapping masses 
indicates good resolution convergence. Halos generally con- 
tain slightly less than the global baryon mass fraction of 
Qb/Qm = 0.13333, indicated as the horizontal dotted line. 
There is a weak trend for smaller halos to contain a smaller 
baryon fraction, as the momentum wind model increases the 
amount of material ejected in these smaller systems. How- 
ever, even the smallest halos that host a resolved galaxy still 
show a modest 30% reduction in the baryon fraction relative 
to the global mean. The stellar fraction shows a more clear 
trend with halo mass, indicating that our momentum wind 
model has more of an effect at suppressing star formation 
than it does at ejecting baryons from halos, as we saw in 
Figure 01 

The middle and bottom panels show the same plot for 
the constant and no wind models, respectively. Unlike the 
momentum wind case, these model shows no trends with 
halo mass in their baryon mass fractions, and a reduced 
trend in their stellar mass fractions. So, for instance, in the 
lowest mass halos, there are more than three times as many 
stars in the constant wind model than in the momentum 
wind case. Again, this can be taken as some measure of 
modeling uncertainty. 



In all cases, halos containing resolved galaxies have 
masses > 5 x 10^ Mq and circular velocities > 50 km/s. The 
circular velocity scale is shown on the top a xis, calculated as 
V200 for a Navarro. Frc nk fc White"! 1^93) profile; the peak 
circular velocity would be higher by 20-30% assuming they 
follow an NFW-like rotation curve. This means that most 
of these galaxies ar e unlikely to be suppresse d significantly 
by photoionization iThoul fc Weinberg Ill99il . regardless of 
whether the photoionization field is accurately modeled or 
not. This is also evident in the baryon fractions; had gas 
infall been suppressed due to photoionization heating, the 
baryon fraction would be significantly reduced in smaller 
systems. We have examined smaller volume (4/i~'^Mpc box 
size) simulations and we do in fact see significant suppres- 
sion of the baryon fraction in halos with Afhaio ^ W^'^Mq, 
which is smaller than any halo with a resolved galaxy in 
these runs. It should be pointed out that these galaxies and 
halos grew from smaller systems that may have been more 
strongly affected by a local radiation field (if it was present) 
in the past, but for particularly the more massive galaxies 
it appears they will have spent most of their growth phase 
in a regime that is broadly unaffected by whether or not a 
local photoionizing radiation field exists. 



4 CLUSTERING AND INHOMOGENOUS 
REIONIZATION 

Galaxy clustering is an important ingredient in the reioniza- 
tion process. Early galaxies form within highly biased peaks 
in the matter distribution, that are highly clustered. This 
may enhance the inhomogeneity of the reionization process, 
and may even result in a characteristic bu bble size for cos- 
molo gical HII regions in the early univers ejAj^rit he fc Loeb I 
200l lFurlanetto. McQuinn. fc Hernguist 120091 . In this sec- 
tion we discuss clustering of our high-redshift galaxies, and 
the implications for local reionization. 

4.1 Evolution of Clustering 

Figure |H| (top panel) shows the two-point correlation func- 
tion from the 16ft~^Mpc momentum wind model computed 
for dark matter {^dd, dashed) and galaxies with stellar mass 
M* « 1O**-^~''M0 {^gg, solid). Thin and thick lines indicate 
2 = 9 and 6, respectively. The redshift evolution of for 
mass-selected galaxy samples is determined by a competi- 
tion between the formation of new galaxies in less biased 
areas and the gravitational growth of matte r clustering. As 
noted at lower redshifts bv I Weinberg et all fe.g. l200'3l . this 
can have the result that ^gg is nearly constant with redshift 
even though the comoving correlation length of the underly- 
ing dark matter increases significantly. This effect can also 
be seen in the bottom panel of Figure IHl w hich gives the 
redshift evolution of the bias b = \/ S.gg/S.dd- As expected, 
the bias of the mass-selected sample drops dramatically be- 
tween 2 = 9 and 2 = 6, with the largest evolution occurring 
at the smallest scales. 

We have fitted power laws of the form ^(r) = (r/ro)"'' 
to the correlation functions for the resolved galaxies and 
dark matter in each of our simulations. Figure |5] gives the 
redshift evolution of ro for the 8, 16, and 32/i~^Mpc vol- 
umes of the momentum wind runs, with approximate mass 
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Figure 8. Top: Two-point correlation function for dark mat- 
ter (daslied) and galaxies with stellar mass M* > 10*^^ Mq 
(solid) from the 16/i~^Mpc momentum wind model. Thin and 
thick curves indicate z = 9 and 6, respectively. Bottom: Bias 
b = i^gg /^dd)^^'^ ^ ^ function of scale in the same model. Galax- 
ies reside in highly biased regions, and the bias is a strong function 
of redshift and scale. 

ranges as indicated in the caption. The dark matter corre- 
lation length is also plotted, multiplied by 10 to increase 
visibility. We have not attempted to measure the correla- 
tion function of the resolved galaxies in the 32fe~^Mpc box 
at 2 = 8 — 9 because this simulation contains too few re- 
solved galaxies at these redshifts. As anticipated by peak- 
bias formalism analyses, early galaxies are highly clustered, 
with a strongly scale-dependent and mass-dependent bias. 
Galaxy correlation lengths evolve only weakly with redshift, 
despite the fact that dark matter correlation length doubles 
from 2 = 9 — > 6. This suggests that the cumulative effects 
of photoionization from collections of galaxies may signifi- 
cantly accelerate the reionization process in local regions. In 
the next section we estimate the strength of this effect. 

4.2 Overlapping Galaxy Stromgren Spheres 

Before reionization completes, the moment at which the ion- 
ized region around a galaxy overlaps with a neighboring 
galaxy's ionized region depends on the rate at which the 
galaxy leaks ionizing photons into its surrounding IGM as 
well as the distance to the next galaxy. In this section we 
estimate the size of the ionized regions and compare this to 
typical distances between galaxies, estimated using a neigh- 
bor search as well as from the correlation function, in order 
to determine whether galaxies' Stromgren spheres overlap 
by 2 = 9. Our calc ulation is similar in spirit to that of 
iBarkanal ll2002l) and iFurlanetto. Hernauist. fc Zaldariagga I 



Figure 9. Top: Comoving correlation length versus redshift for 
dark matter (long-dashed) and resolved galaxy samples from our 
three momentum wind volumes. The stellar mass ranges for galax- 
ies in the volumes in increasing order of size are M* f» lO^'^^^M© 
M* lO^ i^^M©, and iW* ~ W^^'^°Mq. The dark matter corre- 
lation length has been multiplied by a factor of 10 for visibility. 
The error bars are Poisson, and do not include cosmic variance. 
They have been offset slightly for clarity. Despite an increase of 
a factor of two in the dark matter correlation length, the galaxy 
correlation length is roughly constant from 2 = 9^6. Bottom: 
Slope of the correlation function versus redshift for the same sam- 
ples. The slopes are broadly similar to that seen locally. 

i2004l) . with the main difference being that we will use the 
gas distribution taken directly from our simulations rather 
than one based on peak-bias formalism. 

In order to estimate how much of the surrounding IGM 
each galaxy ionizes, we determine how many ionizing pho- 
tons the galaxies had emitted up to 2 = 9 by combining 
the simulated star formation history of each galaxy with a 
fit to the number of i onizin g photons emitted in the stellar 
models of Schacrer] 1)20031 equation 1). These models as- 
sume a Salpeter initial mass function (IMF) between 1-100 
Mq and account for the effects of metallicities down to near 
metal-free populations. Typically, with this prescription we 
find that 1 Mq of stars formed by 2 = 9 ionizes « 10* M© of 
gas. This is more than twice the ionized mass that would be 
expected f or sol ar metallicities (~ 4OOOM0), in accord with 
ISchaererl l)2003fl , which in turn in some what more than tha t 
expected for a standard Salpeter IMF jMadau et alJll999ri . 
From the number of ionizing photons we determine the mass 
of gas that would have been reionized assuming a hydro- 
gen mass fraction X=0.75, escape fractions fesc = 1 and 
0.1, and no recombination (we estimate recombination radii 
later) . Finally, we determine the radius of the sphere around 
each galaxy in the simulation that encloses this mass of gas. 
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Figure 10. Radius of Stromgren sphere versus stellar mass for 
galaxies at z = 9. The y-axes are in comoving (left) and proper 
(right) units. Three-pointed and four-pointed stars denote the size 
of the ionized IGM regions assuming escape fractions /esc = 0.1 
and 1, respectively, while the open squares show the radius of a 
region that can recombine in a Hubble time. Open circles and 
open triangles give different estimates for the overlap radius as 
described in the text. Galaxies less massive than 10^ Mq have 
probably already achieved overlap with a similar galaxy and ion- 
ized their infall regions by z = 9 assuming /esc 0.1, while more 
massive galaxies have probably achieved overlap with less massive 
galaxies. 



including the galaxy itself. This ionized sphere radius repre- 
sents a crude estimate of the extent to which a galaxy can 
ionize its surrounding IGM. 

Fieure lTITl shows. as a function of stellar mass, the size of 
the ionized sphere around each resolved galaxy in the three 
momentum wind runs. Triangles and squares correspond to 
fesc = 1, fesc = 0.1, respectively; as we discuss below, the 
latter is probably more relevant, while the former is included 
to show the sensitivity to this uncertain parameter. These 

1/3 

points closely follow a slope r oc M, . Since the ionizing flux 
is closely tied to the stellar mass, which in turn is reasonably 
closely tied to the halo mass (cf. Figure|7|l this indicates that 
the universe is fairly homogeneous at these radii. In words, 
this plot indicates that galaxies with stellar mass around 



are capable of ionizing spheres around ~ lh~^ Mpc 
(comoving) in size, assuming an escape fraction f^sc 5S 0.1; 
regions ionized by lower-mass galaxies scale in size with the 
galaxy's stellar mass roughly as r oc M^^. 

Do the ionized regions overlap by z = 9? A brute-force 
approach to this question is to determine distance to the 
nearest neighbor for each galaxy in each simulation, take 
the median, and divide by 2. This distance, which we call 
the overlap radius R, is indicated by the large open circles 
in Figure IIUI the error bars indicate the mass ranges that 
the samples cover. A slightly more elegant approach is to 
find the radius 2R of the sphere that is expected to enclose 



one neighbor using the known mean number density no and 
galaxy autocorrelation function: no{$,{r) + l)dV — 1. 
Using this approach, we obtain overlap radii given by the 
open triangles, which are somewhat smaller than that ob- 
tained using our brute-force approach. No estimate is made 
for the largest mass bin because there are too few galaxies at 
2 = 9in the 32/i"^Mpc run to reliably estimate a correlation 
function. 

Comparing the size of the ionized regions to either over- 
lap radius for each sample indicates that galaxies with stel- 
lar masses 1O'^~^M0 have already achieved overlap with a 
similar galaxy by 2; = 9. By contrast, the galaxies in the 
most massive sample are too rare to have achieved overlap 
with a similar galaxy a,t z — 9 despite the large ionized re- 
gions that they can create. By z = S, R drops to 3.5/i~^ 
Mpc for these galaxies so that they do overlap. Neverthe- 
less, since these massive galaxies fall in highly biased regions, 
it is reasonable to expect that by z = 9 they will already 
have achieved overlap with nearby lower-mass galaxies (as 
argued bv e.g. [ Furlanetto. Hernauist. fc Zaldariagga1l2004l : 
IWvithe fc Loeb Il2005^ . This means that their infall region, 
i.e. the region of space where gas is falling onto that galaxy 
as opposed to any other, has likely been ionized by then. 

Of course, the actual physical situation is considerably 
more complicated than the simple model we have presented 
here. First of all, ionizing photons must escape from the 
galaxy as well as the surrounding dense IGM in which the 
ionized hydrogen would be able to recombine quickly. An 
estimate of size of this region is shown as the open squares 
in Figure [101 This is computed as the radius of the sphere 
which encloses the mass of gas (including the galaxy) which, 
had it been ionized very early on, would have been able to 
recombine within a Hubble time. The recombination time 
becomes very small close to the galaxy, so in order for the 
IGM to become ionized, some fraction of the ionizing flux 
from each galaxy must be permitted to escape beyond this 
rad ius. A canonical value fo r this escape fraction is ~ 10% 
(cf. iHaardt fc Madau ll200lD . motivating our above choice 
for this parameter. Our analysis shows that once an ionizing 
photon does manage to escape into the IGM, the recombina- 
tion time becomes comparable to a Hubble time and most of 
the ionized gas can be expected to remain ionized. Recom- 
binations will also reduce the size of the Stromgren sphere 
somewha t, but for these bubble sizes the effect should not 
be large jFurlanetto fc Oh 2005) . 

Further complications that we do not consider in this 
model are the non-sphericity of the matter distribution that 
results in a more complex morphology for ionized regions, 
the fact that the universe is denser at earlier times and 
therefore the recombination region may actually be some- 
what larger, and effects of light-travel time and photon 
redshifting. A full cosmological radiative transfer hydrody- 
namics code is required to model all these effects, which 
is beyond the scope of this paper. We note that a num- 
ber of radiative transfer ca lculations for the early universe 
have already been done (e.g.lGnedin l 2000l:lRazoumov et al. I 

■ "HI 



2002': 'Ciardi. Ferra ra fc White Il200a ISokasian et al. II2004I : 

JQiev et al 2005,) . so we hope to include full modeling of 

the topology of reionization in the future. But for the pur- 
poses of this work, our quasi-analytic modeling is likely to 
be qualitatively accurate. 

In summary, our simulations combined with simplistic 
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Figure 11. Luminosity functions of galaxies in the momentum 
wind model for various observed IR bands at z = 8 (top panel) 
and. 2 = 6 (bottom panel). The luminosity functions are fairly 
steep, and show that galaxies are generally blue, becoming bluer 
at higher redshifts. 
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Figure 12. Top: Evolution from z = 9 ^ 6 of the luminosity 
functions of galaxies in the momentum wind model in the Spitzer 
[3.6/i] band. Bottom: Comparison of [3.6/i] luminosity functions 
in the constant wind and momentum wind models, at z = 8 and 
z = 6. The constant wind model has higher luminosities and a 
steeper faint-end slope. Units of # are number per (/i~^Mpc)^ 
per magnitude. 



modeling suggest that, by z = 9, galaxies more massive than 
IO^Mq are likely to have achieved ionized bubble overlap 
with other galaxies even if the reionization process has not 
yet completed in less biased regions. These galaxies should 
have therefore ionized their infall regions prior to z ~ 9. The 
overlapping Stromgren spheres should contain within them 
a photoionizing b ackground strength that i s comparable to 
that estimated bv iHaardt fc Madau I i200ir) . Hence our as- 
sumption of a uniform photoionizing background, while not 
valid for the bulk of the volume, may be reasonably accurate 
for the regions around the galaxies that we are interested in. 



5 PHOTOMETRIC AND EMISSION LINE 
PROPERTIES OF GALAXIES 

High-redshift galaxies have mainly been discovered and ob- 
served using two techniques: Broad-band colour selection 
and Lya emission line searches. Broad band selection, usu- 
ally utilizing the Lyman continuum edge or even the Lya 
edge at high redshifts, offers the advantage of requiring only 
photometric surveys. However, galaxies are faint in these 
broad bands, particularly reionization-epoch objects which 
require infrared observations that face high sky backgrounds 
from the ground. Also, the precise redshift cannot be deter- 
mined, hampering three-dimensional clustering and environ- 
mental studies. Lya selection has the advantage that sources 
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Figure 13. Comparison between observations o f the rest- 1350 A 
lumin osity function of 2; ~ 6 j-dropout galaxies bv lBouwens et al. I 
J200(i) versus our three wind models. The momentum wind model 
predictions are in good agreement with data, the constant wind 
model predicts about three times too many galaxies, and includ- 
ing no superwind feedback results in more than an order of mag- 
nitude too many galaxies. 

may be very bright in young, dust-poor systems as expected 
in the early universe, they can be found in narrow-band 
filters that offer low sky background, and a fairly precise 
redshift is immediately obtained. However, Lya is a difficult 
line to interpret, as it is easily scattered and extincted, so 
that even a smattering of dust can result in highly attenu- 
ated flux. A full understanding of reionization epoch galaxies 
will likely require pursuing both approaches. In this section 
we present the broad band and Lya emission properties of 
z = 6 — 9 galaxies from our simulations. 

5.1 Photometric Properties 

We examine a series of bands redwards of redshifted Lya 
at z > 6, since bluewards bands are expected to be heav- 
ily attenuated by intervening neutral hydrogen. Figures [TTI 
and II 21 show various luminosity functions (LFs) for the two 
wind models for J [24/i] bands, corresponding roughly to 
rest-frame 1500A to A'-band at z ~ 6 — 9. These LFs have 
been obtained by stitching together the individual LFs from 
our three volumes for galaxies above our resolution limit, 
in a manner similar to that described for mass functions in 
There is generally good agreement in the overlap regions 
resulting in a smooth LF from 24th to 32nd magnitudes, 
although the overlap regions are typically less than a mag- 
nitude. The turnover at faint magnitudes arises from our 64 
star particle mass resolution limit in our smallest volume. 



Figure 1111 shows the z = 8 and z — 6 predicted LFs 
in four infrared broad bands. Galaxies are generally blue, 
with the rest-K LF lower by about two magnitudes at both 
redshifts. Hlence while galaxies possess an old stellar pop- 
ulation, their light is dominated by young stars. Galaxies 
at 2 = 8 are somewhat bluer than at z = 6, as at z = 8 
the rest far-UV LF (observed J) is larger than the rest blue 
LF (observed [3.6/^]), which is larger than the rest optical 
(observed [8/i]) one; by z = 6 the LFs for all these bands 
are similar. Galaxies as bright as 25th magnitude in near-IR 
bands exist at z = 8, but they are quite rare and presum- 
ably highly clustered, requiring large-area surveys to find. 
The steep LFs suggest that modest increases in depth will 
yield substantially larger samples. 

Figure^] (top panel) shows the evolution of the Spitzer 
[3.6/i] band LF from z = 9 — > 6 in the momentum wind 
model. The LF shows a steady increase, by about three- 
quarters of a magnitude per unit redshift, slowing a bit from 
2 = 7— > 6. At 2 = 9 the LF looks like a power law, but by 
2: = 6 it obtains a slight bend, akin to the stellar mass func- 
tions in Figure 13 Similar evolution is seen in other bands. 

The bottom panel shows a comparison of the z = 8 
and z — 6 LFs between the momentum wind and the con- 
stant wind models. The constant wind model, by virtue of 
its greater amount of star formation, has higher luminosity 
functions. The shape is also slightly affected, with the con- 
stant wind model showing a steeper faint-end slope. This is 
a reflection of the stellar fraction varying more with mass 
in the momentum wind model than in the constant wind 
shown in Figure Q the high mass loading factor in 
small galaxies is suppressing the faint-end slope in the mo- 
mentum wind model. This may prove advantageous if the 
trend continues to lower redshifts, since simulated luminos- 
ity functions from constant wind simulations at 2 ~ 4 ap- 
pear to have too steep a faint end compa red to observations 
ilNight et al.1l2005l: iFinlator et aLlbOOSD . 

Observ ations of 2 ^ 6 galaxie s are already rapidly ac- 
cumulating. iBouwens et al. I ll2006ll found over 500 systems 
in Hubble NICMOS and ACS parallel fields that satisfy i- 
dropout Lyman break selection criteria. Of course, some are 
probably interlopers, but they estimate a rather small level 
of contamination (< 8%). Figure IT^ shows a comparison 
of our z — 6 luminosity functions computed at rest-frame 
I350A to their observations. Overall, the agreement in both 
shape and amplitude is quite impressive, particularly for the 
momentum wind model (solid line). Note that no parame- 
ters have been tuned specifically to obtain this agreement. 
The no wind model (dashed line) overpredicts the number of 
galaxies at a all luminosities by well over an order of magni- 
tude, showing that feedback must suppress galaxy formation 
at all masses in order to match observations. iBouwens et al. I 
i2004l) detected a handful of 2 ~ 7 — 8 systems down to 
H < 28, and found that the evolution from 2 ~ 4 is about a 
factor of five in number; we predict (cf. Figure 1121 a factor 
of ~ 3 increase from 2 « 8 — ^ 6, so the observed evolution is 
not inconsistent with our models, though this bears further 
examination down to 2 = 4. 

While this comparison favors the momentum wind 
model at face value, it is worth recognizing that the number 
of high-redshift galaxies is sensitive to the CDM power spec- 
trum on small scales, which is still not precisely constrained. 
For instance, early galaxy formation can be significantly 
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Figure 14. Sample galaxy spectra drawn from the 16h~^yipc 
momentum wind simulation. These generally show a significant 
4000A break and a blue continuum, with some variation due to 
difli^erences in star formation histories. The log of the stellar mass 
of each galaxy is shown on the right. The top one is the most 
massive galaxy in this run at this epoch, and the rest are ran- 
domly selected within various ma ss ranges. The data points are 
observations of a 2 7 galaxy bv lEgami et al.1 l2005ft corrected 
by a factor of 25 for lensing magnification. The data are well fit 
by one of our galaxies with 10^''^ Mq. 



suppressed in warm dark matter (lYoshid a et al. Il2003all or 
running spectral index ( Yoshida ot al. 20 03tJl models. Since 
these ACDM variants rem ain viable at p resent (though in- 
creasingly disfavored: Seliak et al. |l2004l) . we do not claim 
to constrain our feedback models based on z > 6 galaxy 
abundances. Furthermore, the new WMAP3 data favors a 
slightly lower matter power spectrum at these redshift at the 
relevant scales, which would lower our predicted LPs. Hence 
we leave a more quantitative study of this for future work. 
Here we simply note that our models broadly agree with ex- 
isting observations of z ~ 6 LFs, and with improvements in 
cosmological parameter precision such data could play a key 
role in constraining models of early galaxy formation. 

5.2 Sample Galaxy Spectra 

Thanks to unprecedented sensitivities of Spitzer and other 
near-infrared instruments, it has now become possible to 
study the stellar populations of galaxies at z > 6 at some 
level of detail through multi-band photometry spanning into 
the rest-optical. This provides even more stringent tests on 
our simulated galaxy population. Here we examine the rest- 
UV and optical spectra for a fe w of our simulated galaxies . 
These spectra are obtained from lBruzual fc Chariot 1 (l2003l) 



population synthesis models using the Padova 1994 models 
and a Chabrier IMF, where each star particle in a galaxy 
is treated as a single stellar population whose age is com- 
puted from the star particle's time of creation. The single- 
burst stellar population spectra of a given galaxy are then 
summed to obtain that galaxy's spectrum. We interpolate 
the population synthesis models to the stellar metallicity of 
each galaxy. Note that we are not employing a top-heavy 
IMP, since the metallicities of these galaxies fall within the 
range where there are no compelling arguments for devia- 
tions from locally-derived IMPs. 

Figure ITU shows a selected sample of spectra from the 
16/i~^Mpc simulation at 2; = 7. As can be seen, most of these 
systems possess a significant old stellar population as evi- 
denced by a significant 4000A break, and have a contirmum 
that is typically slightly blue through the rest-UV, consis- 
tent with the LFs shown in Figure fT?l There is some varia- 
tion in the continuum shape due to the details of individual 
star formation histories, but overall there is no trend with 
mass, because the birthrates of galaxies are similar at all 
masses (cf. Figure |1J. Indeed, these galaxies look not terri- 
bly dissimilar to present-day star forming galaxies, with the 
exception of their low luminosity (and of course no emission 
lines since these are not present in the Bruzual & Chariot 
models) . 

Recently, lEgami et al. I ^2005) observed a putative z ~ 
7 obj ect found though a cluster lens search iKneib et al. I 
l2004h in Spitzer and Hubble near-lR bands. They showed 
that this system was best modeled as having ~ 5 — 11 x 
10* Mq of stars with a substantial old stellar population. 
Their data, corrected downwards by a factor of 25 for lensing 
magnification, are shown as the data points in Figure 1141 
A galaxy randomly drawn from our 16/i^^Mpc momentum 
wind simulation having a stellar mass of 5 x 10* just so 
happens to be an excellent match to these data. This galaxy 
has a star formation rate of 1.5 Mq/jt, a median stellar 
age of 130 Myr, and one-tenth solar metallicity, which are 
all typical values for a galaxy of this size in our simulations. 
They are also consistent with values inferred from a more 
detailed analysis of this object by Scliacrcr & P cUo . (■2005i) . 
This shows that galaxies like the Egami et al. object are not 
only expected, but are fairly common and typical systems 
at this epoch. 

Another object, a z — 6.56 lensed Lya emitter found 
bv iHu et al. I ll2002D. was follow ed up with Spitzer by 
ICharv. Stern fc Eisenhardt I (|200^ , who found that its spec- 
trum is similar in shape but intrinsically somewhat brighter 
than the Egami et al. object. Not surprisingly, it is possi- 
ble to find simulated galaxies with M IO^Mq that match 
most of the broad-band data for this object as well. How- 
ever, this object does have a peculiarly high IRAC [4.5/i] 
fiux. By fitting a model to the remaining data and assum- 
ing the excess flux is due to unresolved Ha line emission. 
Chary et al. estimate a star formation rate of ~ 14OAf0 /yr. 
This star formation rate would be far above that expected 
for galaxies of this mass (cf. Figure QJ, so if confirmed this 
would present a challenge to our models. However, this ob- 
ject was also analyzed bv {Scha erer & PcUo (200^ prior to 
the Spitzer data becoming available, and their estimated star 
formation rate is 11 — 41M0/yr, much more in line with 
our simulations. Another object, the claimed detection of a 
6 X W^^ Mq system with a photometric redshift of 2 ~ 6.5 
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5.3 Lya Emission 
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Figure 15. Top: The Lyman alpha emitter luminosity function 
at z = 6 in the momentum wind (solid), constant wind (dotted), 
and no wind (dashed) models, assuming a Lya escape fraction of 
2%. The units of <I> are number per unit (h~^Mpc)'^ per log lumi- 
nosity interval, and Li^y^ is in erg/s. Data from Santos 1 2004) are 
shown for comparison. For /esc = 0.02 the mzw model matches 
reasonably well to the data. Bottom: The evolution of the Lya lu- 
minosity function from z = 9 — > 6 in the momentum wind model, 
for an assumed escape fraction of 2%. The evolution is about half 
an order of magnitude in Li^y^ at a given number density over 
these redshifts. 



bv lMobasher et al. I (120051) . would be difficult for our mod- 
els, or indeed any ACDM-based model, to reproduce, as the 
expected space density of such systems is less than one per 
cubic gigaparsec. However, other groups favor a substan- 
tially lower redshift (and lower mass) for this object, so it 
remains to be seen if this galaxy will indeed force revisions 
to our cosmology. 



Lya emission galaxies are now being detected o ut to 

z ~ 6.5 using narrow-band surveys llHn et al. I l2002l : 
iMalhotra fc Rhoadslliooi IStern et al. II2005I) . Higher red- 
shifts move into the infrared, where current space-based ob- 
servatories do not have sufficient collecting areas to detect 
these faint sources, necessitating ground-based observations 
through a bright night sky. Fortunately, there are some rela- 
tively clean, albeit narrow, windows between water and OH 
lines where the sky background is low (see e.g. lBarton et al~l 
|20o3)- A number of narrow-band surveys are underway to 
take advantage of these windows to searc h for Lya emit- 
ters at z ~ 7 - 10 Barton et ar| 2004: Ho rton et al. II2004I : 
IWillis et al.1l2005l: Istark fc Ellis 11200511 . In this section, we 



make predictions for the Lya emitter population from our 
simulations for comparison with such observations. 

As discussed in lBarton et al. I (|200^ , there are large un- 
certainties regarding the detectability of Lya emission from 
star formation. Scattering due to the surrounding IGM in- 
cluding the damping wing may make Lya emission more 
diffuse and difficult to detect. As we showed in H3.4I signif- 
icant metals are expected to be present, so dust extinction 
may be considerable, unless the ISM remains surprisingly 
unmixed. On the other hand, if the initial mass function 
is extremely top-heavy, then perhaps Lya emission can be 
boosted. 

To make predictions for Lya luminosity functions from 
our simulations, we convert the instantaneous star forma- 
tion rate of each galaxy into a Ly a luminos i ty usin g a fit to 
the models listed in Table 4 of ISchaerer I i2003h . In par- 
ticular, we use a solar metallicity Lya emission value of 
2.44 X 10'*^ erg s~^Mg^ yr, and then correct for individual 
galaxy metallicities using equation 1 of ISchaerer I i2003l) . Of 
course, not all Lya photons emitted by young stars will ac- 
tually reach us; there is an escape fraction that is highly 
uncertain. Note that by "escape fraction" here we mean the 
fraction of Lya photons that reach our telescopes, not the 
fraction that escapes from the star forming region or galaxy 
itself that is canonically around 10%. The escape fraction 
depends on the amount of scattering and absorption in the 
ISM and IGM along the line of sight, which is difficult to 
predict directly. 

Here we take the approach of constraining this escape 
fraction empirically by comparing our z = 6 p r edicte d 
Lya luminosity function with data from ISantos I i2004l) . 
They used Keck/LRIS to obtain spectra of objects in high- 
magnification regions of foreground clusters, and found 11 
likely (7 confirmed) Lya emitters with 2; « 4.6 — 5.6. Lensing 
magnification offers the best current technique for detecting 
faint sources, and indeed these objects represent the faintest 
Lya emitters currently known. The downside is that lensing 
results in a considerably smaller volume being probed. 

In Figure IT^ we present simulated luminosity functions 
for Lya emitters from the three wind models (top panel) 
and the evolution from z = 9 ^ 6 in the momentum wind 
model (bottom panel). We choose an escape fraction of 2% 
for Lya photons that would be dete ctable by us, w hich is 
selected to obtain agreement with the|SantosJ (|200^ obser- 
vations (using all 11 likely sources) in the momentum wind 
case, shown as the data points. For the constant wind case, 
an escape fraction of 0.5% is an equivalently good fit, and for 



16 Dave, Finlator, Oppenheimer 



the no wind case the best-fit escape fraction is around 0.2%. 
Note that the escape fraction tunes only the ampHtude of the 
luminosity function, while the slope remains an independent 
prediction of the models (assuming that the escape fraction 
is not luminosity-dependent). Therefore the excellent agree- 
ment from log LLya = 40.5 — > 42 is encouraging. Though the 
Santos et al. data would be more ap propriately compared to 
« ~ 5 predictions, IHu et alH (120051) has demonstrated that 
there is no discernible evolution in the Lya luminosity func- 
tion from 2 ~ 4 — > 6.5, so the comparison to our simulated 
z = 6 sample shou ld be re a sonab ly valid. 

Malhotra & Rhoads I ^00^ and iHu et al. I i2005h 
also argue that the lack of evolution from z ~ 5.7 to 
z « 6.5 would be surprising if the universe was predomi- 
nantly neutral at z ~ 6.5, since the damping wings from 
IGM absorption would significantly reduce the fiux and 
change the line profiles. Even if the bulk of the universe 
is neutral, it is possible that the bright Lya sources seen 
at z — 6.5 may have ionized their loc al surroundings , 
facilitating t he escape of Lya photons iHaiman I 120021 : 
ISantos Il2004^ . In either scenario, it appears Lya emission 
will continue to escape from galaxies and be detectable 
into the reionization epoch. Our crude estimates of the 
local photoionization volume in Figure 1101 supports the 
idea that even out to z ~ 9, IGM damping wing en- 
croachment may not significantly attenuate Lya emission, 
which agrees with conclusions reached by various semi- 
analytic studies llFurl anetto. Hernauist. fc Zal dariagga I 
120041: 1 Wvithe fc Loeb II2005I : iMalhotra fc Rhoads Il2005^ . 

If we take our estimated escape fraction at 2 = 6 and 
assume that it does not evolve with redshift, then there is rel- 
atively little evolution predicted from 2 = 9 ^ 6, as shown 
in the bottom panel of Figure \TE\ This scenario is broadly 
consistent with Hu et al.'s lack of observed evolution to lower 
redshifts. At a given number density, the evolution is around 
a factor of three in luminosity between 2 = 9 — ^ 6, or equiva- 
lently a factor of five in number density at a given luminosity. 
This means that surveys designed to probe 2 « 9 will still 
have plenty of candidate targets, they must "merely" go 3x 
deeper than those at 2 ~ 6. 

One instrument that is being currently developed to 
probe Lya emitters in the reio nization e poch is the Dark 
Ages 2 Lya Explorer, or DAzLE jHorton et al. 2004), set to 
begin operation in 2006. Using an _R = 1000 survey in J-band 
night sky windows on the Very Large Telescope, DAzLE 
will survey 6.8'x6.8' patches down to 2 x 10~^* erg s~^cm~'^ 
in the redshift range 2 — 7.7 — 10. We can make predic- 
tions for the 2 = 7.7 atmospheric window, with a red- 
shift ran ge of Az = 0.05 ove r which the sky is relatively 
dark (see iBarton et al. II2OO4I for a diagram of the J-band 
night sky spectrum). Here, their flux limit translates to 
lO*^ '' erg/s, and the window's volume for the DAzLE field 
is 2500 {h~^Mpc)^ . Our predicted source density above this 
fiux limit is approximately 0.002 (/i~^Mpc)"^, so DAzLE 
should expect to see around 5 sources per field within this 
night sky window. The numbers are similar for the 2 = 8.2 
window, which is a bit wider but further away, and the 
2 = 8.8 window which has Az ~ 0.1 but half as many 
predicted sources than at z — 7.7. These numbers are not 
large, but with a few adjacent fields it may be sufficient to 
do preliminary clustering studies to possibly constrain the 
patchiness of reionization. If no objects are seen, this would 



suggest that the escape fraction is dropping substantially to 
higher redshifts, whereas an abundance of objects would sig- 
nify a dramatic change in the IMF or dust content of 2 > 6 
systems. More accurate predictions for detailed comparisons 
will require the development o f Lya radiative tr ansfer codes 
including dust extinction (e.g. iTasitsiomi Il2006^ . 

Another key design issue is the expected line widths of 
Lya emission. Since halos at this epoch have circular veloc- 
ities typically less than 150 kms^^ (cf. Figure|7|l, the naive 
expectation is that emission line widths will not be signifi- 
cantly larger than this. It may be significantly smaller if the 
emission region (or last scattering surface) is not sampling 
the full p otential well, as seen e.g. for Lyman break galaxies 
at 2 ~ 3 jPettini et al. 1120011) . Conversely, it could be larger 
if outflows stimulate Lya emission, or if radiative transfer 
effects smear out the emission region beyond the galaxy's 
host halo, though these scenarios seem relatively improba- 
ble. Hence it is likely (but not guaranteed) that an instru- 
ment with R — 1000 such as DAzLE would not resolve out 
Lya emission lines, which makes internal kinematic studies 
difficult but is preferable for source detection. 

In summary, the Lya emission properties of 2 > 6 galax- 
ies are an interesting way to constrain the topology and ef- 
flciency of reionization. Our simulations suggest that Lya 
emitters will be observable out to 2 ~ 9 with current tech- 
nology, assuming that galaxies sufficiently ionize their sur- 
roundings so that the escape fraction is similar to that at 
2 < 6. Detection of these objects appears feasible in the 
near future, which should provide crucial insights into the 
reionization epoch. 

5.4 Next Generation Instruments 

A key science goal of JWST is to image the objects respon- 
sible for reionizing the universe. The NIRCam team has a 
gal axy formation prog ram designed to image 2 > 7 galajc- 
les llRieke et al. 112003^ . The survey will have 50ks exposures 
in six near-IR fllters, and lOOks in their 4.4/im fllter. The 
expected depths achieved will be roughly 1 njy, or 31.4 AB 
magnitudes, across these bands. They then plan to do follow- 
up emission line studies on selected high-2 objects. 

Our simulated luminosity functions are coincident ally 
resolved down to about 1 nJy, so we can make general pre- 
dictions for the types of objects that will be observed by 
such a survey. For starters, the predicted number density 
will be large: At 2 « 8, each 2.3'x2.3' fleld will contain over 
400 galaxies per unit redshift (cf. Figure ITTl . which should 
be stra ightforwardly identiflable as Bl-band dropouts (see 
iRieke et al. I2003l for filter set). Note that while large, these 
numbers are still s maller than previou s estimates: The semi- 
analytic model of iBarkana fc Loeb I (|^0Q) predicts many 
thousands of galaxies down to 1 nJy. The main difference is 
our suppression of star formation through superwinds, which 
appears to be necessary in order to match observed 2 ~ 6 
luminosity functions (cf. Figure [T^ . We also predict galax- 
ies will be strongly clustered, but as each NIRCam field is 
over 4/i~^Mpc (comoving) on a side which is longer than 
the correlation length (cf. Figure the field-to-field varia- 
tions should not be too large. The physical properties of the 
observed galaxies are likely to be what we have described 
so far: Mildly enriched systems living in halos of masses 
> with the smallest stellar masses observed be- 
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ing a fewxlO'^MQ. It is not expected that even the faintest 
of these objects will push into the regime of predominantly 
metal- free star formation with little or no dust, though such 
objects may enter into the survey if their stars are excep- 
tionally bright. In short, these are not the primeval galaxies 
you're looking for. 

Turning to Lya emitters, the expected faintness of these 
systems requires larger collecting areas, and fortunately 
work is already underway on the 22m Giant Magellan Tele- 
scope jhttp : //www ■ gmt o . org| | , among others. Such a tele- 
scope with an _R = 3000 near-IR tunable narrow-band capa- 
bility would_2rob£an_order of magnitude fainter than DA- 
zLE llBarton et al. lEo04l) resulting in ~ 20x more sources, 
equivalent in space density to ~ 30th magnitude objects. 
Such spectral resolution may enable kinematic studies as 
well. Such capabilities promise to move Lya emitter searches 
from a handful of discovery objects into the realm of large 
statistical samples for detailed studies. 



6 SUMMARY 

We have presented results for the physical and observable 
properties of reionization epoch galaxies from cosmological 
hydrodynamic simulations. Our simulations naturally pro- 
duce an early epoch of star formation that yields significant 
numbers of reionization-epoch galaxies, in broad agreement 
with available observations. We find that: 

• Systems with 10** Mq in stars are common by 2 ~ 9, and 
by 2 = 6 systems with 10^*^ Mq have a similar space density 
to luminous red galaxies today. The mass functions have a 
steep faint end initially, but become somewhat shallower in 
time, particularly in our favored wind feedback model. 

• Star formation rates are well correlated with stellar 
masses, such that the stellar birthrates show almost no trend 
with mass at these epochs, modulo the impact of stochastic- 
ity in star formation caused by merger-related processes. At 
2 = 8, the stellar mass doubling time for galaxies is around 
0.2 Gyr, while at 2 = 6 it is ~ 0.4 Gyr. Star formation rates 
at 2 ~ 9 ^ 6 range from hundredths to tens of solar masses 
per year. 

• The large stellar masses imply significant metal enrich- 
ment in these systems. While they are generally subsolar 
and can be comparable to the lowest metallicity dwarfs to- 
day such as I Zw 18, their metallicities are still well above 
the putative threshold for ushering in a more normal IMF 
as compared to a top-heavy IMF proposed for Population 
III stars. We predict few if any near metal-free stars forming 
in these systems. 

• The halos containing these galaxies have masses well 
above 10^ Mq , meaning that their virial temperatures exceed 
the primordial cooling limit oi W^K, and hence atomic lines 
are the dominant coolant. Their baryon content is somewhat 
below the global average due to feedback, and different feed- 
back recipes can produce factors of three differences in the 
stellar content, depending on mass. Their circular velocities 
are > 50 kms"'^, implying that these galaxies are unlikely 
to be significantly affected by photoionization. 

• Galaxies at this epoch are highly biased and highly clus- 
tered, with a strongly mass-dependent clustering strength 
and scale-dependent bias. The rapid growth of matter clus- 
tering is offset by a rapidly dropping bias, producing a rel- 



atively constant comoving correlation length. Galaxies are 
sufficiently clustered that they are likely to lie in regions 
that have been locally reionized at a much earlier epoch 
than the time of final cosmological H II region overlap. 

• The broad-band properties of these galaxies show that 
they are generally fairly blue systems, getting bluer to higher 
redshifts. Hence the optimal detection band will be the one 
that is just redwards of Lya. There is evolution of about 
0.5 — 0.75 magnitudes per unit redshift in near-IR bands. 
Predicted number densities vary by ~ x 3 depending on the 
wind model, but the evolution is generally similar. 

• Predictions of Lya emitters are less robust because of 
large uncertainties in Lya transfer out of these systems. As- 
suming a Lya emission detectability similar to that observed 
at 2 ~ 5, then surveys with expected sensitivities of fu- 
ture 20-30m class telescopes will find an equivalent number 
density of 2 = 9 objects as near-infrared surveys down to 
~ 30st magnitude AB. Hence Lya narrow-band surveys can 
be a powerful tool for obtaining large samples of reionization 
epoch objects. 

• Our simulations show good agreement with observa- 
tions of 2; > 6 objects. Individual 2 ~ 7 galaxy broad band 
spectra have been observed (e.g. lEgami et al. I EoOSh that 
are well-matched by simulated galaxies, showing the requi- 
site continuum slope and 4000A break. Ou r 2 « 6 luminos- 
ity function agrees well with observations bv lBouwens et al~l 
^2006), particularly for our momentum wind model (modulo 
uncertainties regarding the exact amplitude of the matter 
power spectrum). The Lya luminosity f unction at 2 » 6 is 
in good agreement with observations bv lSantosI i2004ft . for 
a reasonable Lya escape fraction of 2%. We look forward 
to comparisons with future data to provide more stringent 
constraints for our models. 

The outlook for detecting large numbers of reioniza- 
tion epoch objects is promising. Current Hubble data can 
achie ve depths of 28.5 in near-IR bands llBouwens et al. | 
I2OO6'), albeit in small fields. In addition to Webb, several 
future space-based telescopes have been proposed that will 
do large-area near-IR surveys, including some versions of 
JDEM (Joint Dark Energy Mission), as well as the proposed 
Galaxy Evolution Probe (R. Thompson, PI). Though these 
missions remain firmly in the distant and uncertain future, 
the scientific richness of deep wide near-IR surveys ensures 
that they are likely to be accomplished at some point. On 
the Lya emitter side, a number of operational instruments 
should begin detecting a handful of systems in the imme- 
diate future, but surveys in near-IR night sky windows on 
future 20-30m class telescopes should bring these studies to 
maturity. The timescales for such facilities are of the same 
order as the space-based infrared ones, so it will be an inter- 
esting race to see which technique matures first. The future 
looks bright for lifting the veil of reionization and probing 
the final frontier of observational cosmology. 
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